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MONTE CARLO ANALYSIS OF RAREFIED-GAS DIFFUSION 
INCLUDING VARIANCE REDUCTION USING THE 
THEORY OF MARKOV RANDOM WALKS 
by Morris Perlm utter 
Lewis Research Center 

SUMMARY 

Molecular diffusion through a rarefied gas is analyzed by using the theory of Markov 
random walks. The Markov walk is simulated on the computer by using random numbers 
to find the new states from the appropriate transition probabilities. As the sample mol- 
ecule during its random walk passes a scoring position, which is a location at which the 
macroscopic diffusing flow variables such as molecular flux and molecular density are 
desired, an appropriate payoff is scored. The payoff is a function of the sample mol- 
ecule velocity. For example, in obtaining the molecular flux across a scoring position, 
the random walk payoff is the net number of times the scoring position has been crossed 
in the positive direction. Similarly, when the molecular density is required, the payoff 
is the sum of the inverse velocity of the sample molecule passing the scoring position. 
The macroscopic diffusing flow variables are then found from the expected payoff of the 
random walks. The confidence limit interval on the numerical results obtained for the 
flow variables can be reduced by increasing the number of random walks simulated or 
by reducing the variance of the random walk payoff. By reducing the variance, less 
random walks are needed for the same confidence limit interval; thus less computer 
time is needed for the calculation. Biasing functions are found which change the transi- 
tion probabilities such that the variance is reduced. Both a zero-variance and a 
minimum -variance biasing function are found. The zero-variance biasing function is 
approximated in the numerical calculations which were carried out to illustrate the vari- 
ance reduction procedure. The resulting saving in computer running time due to the 
biasing procedure is shown. 



INTRODUCTION 


The present analysis treats the rarefied diffusion of molecules of one species, evap- 
orated or emitted from one plate, through another nondiffusing molecular species en- 
closed between parallel walls (see fig. 1). The problem has been treated by straight- 
forward analog Monte Carlo without variance reduction in reference 1. The present 
analysis treats the diffusing molecule as a Markov random walk, and the local macro- 
scopic properties are found as the expected value of a random variable, the random walk 
payoff. By biasing the transition probabilities and changing the collision payoffs, we can 
retain the expected Markov walk payoff but reduce its variance so that the Monte Carlo 
result will have a much smaller error. 

Monte Carlo methods have been extensively used in nuclear reactor radiation shield- 
ing problems, and sophisticated methods of analysis have been developed to minimize 
error and reduce computing time (refs. 2 and 3). These Monte Carlo techniques and 
methods of variance reduction would apply to any Markov random walk problem. Similar 
techniques to those used in the present analysis can be applied to thermal radiation and 
other rarefied-gas problems. 

In the present analysis a zero-variance and a minimum -variance biasing function 
are found. These biasing functions apply only for positive event payoffs. However, by 
appropriately defining the scoring procedure, the event payoffs can be made positive. 

The minimum -variance biasing function does not give a zero variance, as the zero- 
variance biasing function does, and so must give a relative minimum variance. How- 
ever, in the Monte Carlo calculating procedure, when only a finite number of events can 
be calculated, the bias functions which would be optimum in the actual calculation have 
not been determined. 

The zero-variance bias function is approximated so as to simplify the Monte Carlo 
calculations. These approximations which reduce the variance are the ’’next event" 
calculation, in which we score the probable event payoff after each event whether or not 
the sample molecule actually reaches the scoring position. Also used is "birth biasing, " 
in which the expected initial or birth event payoff, which can readily be calculated, is 
scored for the initial event irrespective of the initial state of the molecule. Also used 
is "survival biasing, " which prevents the sample molecule from being absorbed after 
only a few collisions, since then it would not contribute to the payoffs of the later colli- 
sion. Survival biasing can be achieved by biasing the transition probabilities so that the 
sample molecule is not permitted to reach the absorbing wall. The large variation in 
sample molecular velocity contributes a large amount to the payoff variance; contribution 
can be avoided by "velocity biasing. " Here we score the expected value of the event 
payoff averaged over all velocities due to the collision, independent of what the actual 
value of the velocity of the sample molecule is. 
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Since in survival biasing the sample history does not end, "Russian roulette" is 
used. When the weighting function is sufficiently small, the random walk is either 
ended randomly with a probability p, or the weight is increased by 1/p. Numerical 
results are given to illustrate the amount of variance reduction obtained by these bias- 
ing procedures. 


ANALYSIS 

As shown in figure 1 the Markov random walk of the diffusing molecule can be rep- 
resented by the sequence of states {Xq, Xj, Xg, . . .}. The symbol X R denotes the 
state of the molecule and refers to a point in position and velocity space {Y n , V R } taken 
on by the diffusing molecule immediately after the n^ collision. The probability cor- 
responding to this random walk is given by 

HtX^Xj, . . .)dX Q dXj . . . = E^Xq^xJXq^X^Xj), . . . dX 0 dXj . . . (1) 

The initial probability Eq(Xq)cK 0 is the probability of the molecule originating in 
cIXq at Xq, where X Q refers to the initial position and velocity of the molecule, which 
in this case is given by the initial position Yq and initial velocity Vq. The transition 
probability K(X n+1 |X n )dX n+1 is the probability of a particle that is at X n immediately 
after the n th collision reaching dX n+1 at X n+1 immediately after the next collision. 
The normalizing conditions that apply are 


= i 

y K (Xn|X n .i)dXn = 1 


( 2 ) 


where the integration is over all X. 

Immediately after each event, n = 0, 1, 2, . . . , where by event we mean the estab- 
lishment of a value X„, we record an "event payoff" P„. This event payoff designates 
the contribution to the desired answer by the molecule immediately after the n 1 event 
has occurred and the value of X R has been established. The event payoff can be a func- 
tion of the entire past history (Xq, Xj, . . ., X n }. 

The random walk payoff for the random walk process starting from the initial 
n = 0 event is then given by 
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( 3 ) 


%= E P n 

n=0 

The expected payoff of the random walk is given by 


= * ■= C C W 8 <*0 

JX 0 • /X n 


dX, 


■/■■■/ (t, P nW X 0> K < X ll X 0> • • • «o dX l 
• / X Q -/X n \n=0 / 


or 


A — Aq + Aj + A2 + • . . 


(4) 


(5) 


where is the expected event payoff 

X k = f' f ( P k )E 0 (X 0 )K(X li X 0 ) • • • K(X kl X k-l )dX 0 dXi • • • dX k (6) 

We can write the transition probability to include l events as follows 


k(0 < X „^ l X „> - / / K < X „ + ll X n> • ■ ■ K <V( • • • ‘“n*-! 


(7) 


If we define K^(XlX') = 6(X - X'), we can obtain Green's function 


G(X|X') = Yj K^(X|X') (8) 

1=0 


Given that an event occurs at X', G(x|X')dX is the expected number of events experi- 
enced by the particle that occurs in dX near X. If we assume the event payoff is only 
a function of its state and not a function of its previous history, the event index of P 
can be dropped so that P n (X) = P(X). Then equation (4) becomes 
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( 9 ) 



p ( x n ) E 0 ( x o) K(n) ( X nl X o) dX 0 


= P(X)E 0 (X’)G(XjX’)dX' dX 


We can define the probability density for state X after n events as 

E n (X) -/ E 0 (X')K^ n \x IX'JdX' 

Letting E(X)dX be the expected number of times a molecule is in state dX at X 
averaged over all initial states, we can write 


°o - /» 00 

E(X) = Yj E n (X) = / E 0 (X’)G(x|X')dX» = / Eq(X’) ^ K^(x|x')dX’ (10) 

n=0 J ^ n=0 


We can now write the expected payoff as 

X=f P(X)E(X)dX (11) 

From equation (10) we can write 

E(X) = E q (X) + f E 0 (X 0 )K(XjX 0 )dX 0 + E 0 (X Q )K(X’ |X Q )K(X |X’)dX* dX Q 

+ JJ E 0 (X 0 )K (2) (X'|X 0 )K(X|X»)dX’ dX Q + . . . 

= E q (X) + fj E 0 (X 0 )K(X|X')[6(X’ - X Q ) + K(X’_|X 0 ) + K (2) (X' |X Q ) + . . .] dX’ dX Q 

= E 0 (X) + J' E(X’)K(x|X’)dX’ (12) 

This result is the well-known Boltzmann transport equation. 
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Monte Carlo Process 


In the Monte Carlo process, we first randomly choose an initial sample molecule 
position and velocity from Eq(Xq). We then score P(Xq). We then randomly choose a 
new velocity and position after collision, X^ from K(XjJXq). We then score P(Xj). 

The particle history is continued in this manner until the history is completed (i. e. , 
the molecule is incident on an absorbing wall). The random walk payoff is then given by 
t]q = P(Xq) + P(Xj) + . . . . This process is repeated for N samples to give the sample 
expected payoff, 

_ N 

Hd-iZ V 0, i w = x (13) 

N i=l 

This gives the desired macroscopic result. By the central limit theorem, the probabil- 
ity & that the difference between the Monte Carlo and the theoretical result is less than 
some error e, called a confidence limit, is given by 

/•e/[a(rj 0 )/VN] ^ 

- e M < e ] / e_t /2 dt ( 14 > 

o 

where a (tjq) is the variance of the distribution of t?q. For 95 percent confidence 
= 0. 95), the confidence limit is given by 


e 



(15) 


2 2 
The variance a (?]q) can be approximated by the sample variance S (^q), which can be 

obtained in the Monte Carlo calculation from 


N 

s!( ’o> = 7-: E Si'V 2 < 16 

w ‘ L i=l 

n 

We can see from equation (15) that the smaller the variance cr (t?q), the smaller the 
error in the Monte Carlo calculation. 

The average computer time to run a sample history a is given by a - T/N, where 
T is the total computer time needed to run N sample histories. We can ratio the con- 
fidence limits for the biased process eg to that for the unbiased or analog process e^ 
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and obtain 


£ B _ g B^0^ 
e A ff AtyP 

where the B subscript refers to the biased process and the A subscript to the unbiased 
or analog process. 

We can compare the computing time needed to complete the random walks for the 
same confidence limits, e B = and obtain 

^B _ ^B CT B^V 

T A 

Then the percent change in the computer running time for the biased case compared to 
the analog case for the same confidence limits can be written as 


Conditional Expected Payoff 

We can find the variance of the random walk payoff a (?]q) as follows. The proba- 
bility density of a random walk from some point { x p X j + i> . . . is given by 
H(X i+1 , X i+ 2 , . . . Ix^. The payoff of this random walk is given by 
7 i ^ = P(Xj) + P( X j + j) + . . . . We can then write the conditional expected payoff for a 
random walk starting at X i as W^(X.) and this is given by 

WitXj) = /mxp = P(xp + f P(X i+1 )K(X. +1 | Xi )dX i+1 

+ /p(X i+2 ) K(2) ( X i + 2|Xi)dXi + 2 + . . . 

PfX^jK^^X^JX^dX^ = / P(X)G(X (X i )dX (18) 




x 100 = 1 - ~ CTb ^ 0 ' x 100 
° A a A^0^ 
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We can see that, for the case where the payoff is independent of the past history, the 
conditional expected payoff will be a function only of its state and not of the number of 
collisions, so that the collision index number i in W^X.) can be omitted. We can see 
that 

w (X i+ i) = P(X i+1 ) + f P(X i+2 )K(X 1+2 |X 1+1 )dX i+2 + f P(X lt3 )K< 2 >(X i+3 |x i+1 )dX i+3+ . . . 

(19) 

By multiplying equation (19) by K(Xj jjx.) and integrating, we can combine the result 
with equation (18) and obtain 

W(X t ) = P(X t ) + / W(X i+1 )K(X i+1 |X i )dX. +1 (20) 

This can be written for the present case as 

W(X») = P(X’)+ J W(X)K(X |X’)dX (21) 

A simpler method of obtaining equation (21) is to write = P(xp + 7 ? i+1 - We then find 
the conditional expected value of 7^ to be 

W(Xj) - P(X,) + f |x i+1 )K(X. +1 |x.)dx i+1 

= P(X i )+y'w(X i+1 )K(X. +1 |X 1 )dX. +1 (22) 

which is identical to equation (21). We can see that the expected payoff can be obtained 
from the conditional expected payoff as follows 

X = *(tj 0 ) = J «f(r ?0 |X 0 )E 0 (X 0 )dX 0 = f W(X 0 )E 0 (X 0 )dX 0 (23) 

or from equation (18) 


x-/e 0 (X 0 ) / P(X)G(X|X 0 )dX dX 0 
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Conditional Expected Square Payoff 


The conditional expected square payoff Q^(xp for a particle leaving can be ob- 
tained as follows: The square of the payoff for a random walk starting at X^ can be 
written as 


- [P^) + o i+1 ] 2 

Then 

Q i (X i )^(, 1 2 |X i ) = P 2 (X i ) + 2P(X i ) f ^ Ix^lK^JX^dX^ 

■ + / ^(’lf + ll X i + l) K ( X i + ll X l) <K i + l < 24 > 

This can be rewritten as 

<W = p 2 < x i> + 2 p < x t> / w ( x it i> K ( x i + il x i> dx i + i + / Qi + i‘ x i + i) K ( x i + il x i> dx i + i 

(25) 


Combining equation (25) with equation (22) yields 

Q^) = 2P(X.)W(X.) - P 2 (X.) + J' Qi+^X^^^X^JX.JdX^j (26) 

Again we see that, for the case where the payoff is independent of its past history, the 
conditional expected square payoff will be a function only of its state and not of the num- 
ber of collisions, so that the collision index number i in Q^(X^) can be omitted. Then 
equation (26) can be written as 

Q(X') = 2P(X’)W(X’) - P 2 (X’) + y Q(X)K(X|X')dX (27) 

As in the relation between equations (19) and (21) we can write equation (27) as 


QO 

« x i> = I - p 2 ^] wjw™ 

n=0 


(28) 


which can also be written as 
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Q(X') = J |^2P(X)W(X) - P 2 (X)] G(X |X')dX 
We can now write the expected squared payoff as 

4o) = f l x o) E o( x o> dx o = f '«V«o(V«o 

= yy*[2P(X)W(X) - P 2 (X)] G(X IX^E^X^dX* dX 
From equations (10) and (11) the variance can now be written as 

= 4o) ■ [«<*0>] 2 = / Q(X 0 )E 0 (X 0 )dX 0 - X 2 

I 

= y* [2P(X)W(X) - P 2 (X)] E(X)dX - [y* P(X)E(X)dx] 


Biasing the Random Walk 

We wish to reduce the variance of the random walk payoff so that the error in the 
Monte Carlo calculation as given by equation (16) will be smaller. We can achieve this 
by biasing the random walk process and at the same time changing the event payoff so as 
to keep the expected payoff for the new process unchanged but at the same time have a 
reduced variance. We can bias the probabilities so that the probability density for the 
biased random walk is 

ff^X^X^Xg, . . .) =e;(X 0 )kJ(X 1 |X 0 )kJ(X 2 |X 1 ) . . . (32) 

where the asterisks refer to the biased results. The subscripts on the biased transition 
probabilities refer to the number of collisions since the initial event and show there 
is a dependency of the biased transition probability on the event number. 

We can distort the event payoff for the random walk as follows: 

t?q = Pj + Pj + V* 2 + ■ ■ ■ (33) 

where 
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P( x 0> 

L 0 (X 0 ) 


A 



P(x x ) 

L 0 (X 0 )L 1 (X 0 ,X 1 ) 





P09 

^(X^L^X^XpLgCXpXg) 


J 


and where the biasing ratios are given by 


( 34 ) 


SVn-l)' 


KfXjXn.j) 


n^O 


L 0 (X 0 ) 


Eogp) 

K 0 (X 0 ) 


n = 0 


The normalizing conditions also must apply 

y E ; dx 0 = y e 0 (x 0 )l 0 dx 0 = 1 

y ^(*»l x n-i)®n -f K < X J X n-l> L n ^n = 1 " = i, 2, . . . 

These conditions are equivalent to the condition 

f • f e 0 (X 0 )L 0 K(X 1 |X 0 )L 1 . . . K(X n |X n . 1 )L kl dX 0 dX 1 . . . dX n = 1 

n = 0, 1, 2, . . . 


(35) 


(36) 


(37) 

(38) 


(39) 


The biasing ratios are subject to the condition that L n (X n X ) -> 0 and Lq(Xq) > 0 
since the biased probabilities must be positive. We can see the equivalency between 
equations (38) and (39) since we can write 

y f E S( x 0 )Kj(Xj I x 0 ) • • ■ K* n (X n |X n _ 1 )dX 0 dXj . . . dX n 
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and 


f f e;(x 0 ) k * 1 (x 1 | x 0 ) . . . < tl (x ntl l x „)® 0 ax, . . . dx n+1 

For these two equations to be true, equation (38) must be true. Then by integration of 
equation (39) over X^, Xg, . . . , we also prove that equation (37) must be true. Thus, 
equation (39) ensures that equations (37) and (38) are both satisfied. 

We can see that the expected value for the biased walk payoff is given by 

X* .= f PjE5dX 0+ fj PlE^KjlXjlX^dX 


+ fJJ P 2 E 0 K l( X ll X 0> K 2< X 2l X l> dX 0 dXj dX 2 + 
1 + ^(^) (E o L o> K ( x il x o> L i «<) ®i + • • • 


= \ 


(40) 


So that the expected biased walk payoff equals the unbiased value. 

The expected biased walk payoff beginning at X^ is given as in equation (18) by 


"'i “ P 1 * J" P !+l K l+l (X i+ll X i) dX i+l 


+ ff p i + 2 K r + 2( x i + 2l x i + l) K i* + l( x i + ll x i) dX i + l <^+2 




i+l ^i+l + • • • 


W i (L i ) 


L 0 L 1 • • • L. 


(41) 


Similarly to equation (24), we can write 
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( 42 ) 


Q* = P* 2 ( Xi ) f 2P*( Xi ) / W* +1 K* +1 (X 1+1 | Xi )dX ltl + ^ +1 K* +1 (X ltl | Xi )dX i+1 


From equations (41) and (22) 


*, X W(X.) P i^ X i) f W ( X i+ l) 

W*^) = L_ = i — - — + / K(X X X )L, . dX 

L 0 --- L ! V" L i y V-- 1^ 1+ 1 1 1+1 1+1 


or 


> (43) 


W*( Xi ) = P*( Xl ) + yw? +1 (X i+1 )K 1+1 (X 1+1 |X 1 )d Xl+1 


Then equation (42) can be written as 


% - [2P*( Xl )W*( Xi ) - P* 2 ( Xl )] + y ^ tl K* tl (X i+1 | Xi )dX i+1 (44) 

Then substituting for 0? ^ in equation (44) by letting the subscript i be i + 1, we ob- 
tain 

Q? = [2P*( Xl )W*( Xi ) - P* 2 (X,)] *f [2Pt +1 (X itl )W* +1 (X i+ ,) - P* + 2 !(X i+1 )] 


XK J + l< X i + xl X i> dX i + l 


(45) 


Continuing this process gives 

$ = ( 2 P*W* - P* 2 ) */ (2P* +1 W* +1 - p *+ 1 )K* +1 (X i+1 | Xi )dX 


i+1 


(L 0 L i • • • L 


*//KA • P U2) K t t l< X i*ll Jt i> K ‘( X i*2l X i *lX K t + l “,*2* 

(W - pf) t /fe± istL^teu!^ 

J L i+1 


■If 


2P i + 2 W i + 2 ' P i + 2 


L i+l L i+2 


K(Xi +2 |Xi + i)K(Xi + i|X i )cJXi + i dX i + 2 


(46) 


We can now write the expected square biased payoff as 


13 




This is th? function we wish to minimize subject to conditions in equation (39). 

A u: 1 simplification is to write the biasing ratio in terms of a biasing function 

I (ref. rs follows: 


L (X ,X ,) = 
n v n’ n-1 7 


W 

W X n-l> 


(48) 


The biased random walk then is 

H*(X 0 ,X 1 , . . .)=[Ej(X Q )][iq(X 1 |x 0 )][l^(X 2 |X 1 )] . . . 


r 

* , , M x i)’ 



[■wvvl 

KX, X Q -i — L 

1 0 W 


k(x 2 x,)-i — L 
2 1 ^(xp 


The normalizing conditions given by equations (37) and (38) are 


(49) 


f 


(50) 


/’ 


K < X „ X „-l> 


I (X ) 
n' n' 




dX = 1 n = 1, 2, 3 


Or equivalently as was shown previously 

/"W = 1 " = 0 . 1. 2 , • • 


(51) 


(52) 
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The payoff for the biased random walk now becomes 


* 

Vq 


_ T>* 
" *0 


+ Pi 


P(x 0 ) P(x x ) P(X 2 ) 

w + w + w 


(53) 


Notice that the biased payoffs are now only functions of the immediate state and do 
not depend on past history. However, they do depend on the number of previous colli- 
sions since the biasing function I R depends on the collision index number n. 

As in equation (18) we can write the conditional expected biased payoff as 


W*( Xi ) = ,(„* |X t ) - P‘(X,) + f P* +1 (X i+I )K* +1 (X i+1 |Xj)dX 1+ j 


f- 

w J 1 


/ P L2< X i + 2> I 4+2< X i + 2 l X i> dX i + 2 


IS±iL K ,X. JxA^dX. 

(V \ 1+1 1 T /V \ 1+1 

i-t-r i+P 


W 


■/ 


P(X i+2 ) ( 2 ) 

+ # — K ^ (X- 


W X i+2> 


i+2 


X i) 


L 


i+2 


w 


W(x.) 

w 


(54) 


Similarly, the conditional expected square biased payoff for the random walk leaving X i 
is given by equation (46) as 


Q*(Xj) = #(^ 2 |x ( ) - 



(n)/ 


w 


2 “ ( x i+n l X i) 

)W(X. ) - P^(X. nI 1+n 1 

l+n' ' i+m ' 1 


, K v 

: i + n>] — 


(X. ) 

i+n v i+n' 


dX. 

l+n 


(55) 


The expected square biased payoff now becomes 


15 



* (»o 2 ) = 


oo 

Jf W Li [ 2PlX " )W(X " 

^ , n=0 v. 


- P 


, . K <n, (X |X„) 

! ( x n )] — ck 


W 


VdX, 



w 

w 




(56) 


where variance for the biased process is given by 

« 2 (,5 2 ) - ^(»5 2 ) - £ 2 (d 5> = ,(,*o 2 ) - * 2 (57) 


Zero-Variance Bias Function 


We wish to choose a biasing function that will minimize this variance. An obvious 
zero-variance biasing function is given by ' 


, P < X n> 

I* n (X n ) = —2- 


if P > 0 


n 


(58) 


where the dagger superscript is used to designate this particular zero-variance biasing 
case. We can see that I* satisfies the normalizing conditions given by equation (52). 

1 n ir 

With I* as the zero biasing function, the event payoff after the n 11 event is given by a 
constant A , regardless of the state X n of the molecule. Thus, for any random walk, 
irrespective of the states involved in the random walk, the individual random walk payoff 
is always given by 


7]q = Aq + A^ + A 2 + . . . A n + ... — A (59) 

where there are an infinite number of terms. Since the random walk payoff is always 

2 t 

equal to the expected random walk payoff, the variance a ( tiq ) must be zero. This can 
also be seen from the previous equations by substituting I^(X n ) from equation (58) 
into (56) to obtain 
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\f [ 2W < X n> - P (V]W®„ 

X „ £ [ P < X n> * 2 £ P < X „ + l> K < X n + l IV^n+l + ■ • •] E „< X „> dX „ 
X n X n+1 


4o 2 )= E X n[ X n + 2 <Vl + X n + 2 + ' • •>] 

n=0 

°° 00 00 / 00 \ 2 

= E X n + 2 E X » E Vi - fE X ») = ^ (60) 

n=0 n=0 i=l \n=0 / 

and from equation (57), a ^(r}Q^j = 0. This result states the obvious fact that a biasing 
function that minimizes the variation in the event payoff to zero will reduce the variance 
of the Monte Carlo result to zero. 

If the event payoff is negative, the zero-variance biasing function would not apply 
since, from equation (58), P > 0. However, it is possible to define a payoff function so 
that it is always positive. For instance, in the scoring of the molecular flux, the sample 
molecules passing in the positive direction would score a positive quantity, while the 
sample molecules in the negative direction would score a negative quantity. Then, since 
part of the payoff is negative, the zero biasing function would not apply. However, we 
can score the sample molecules in the positive direction and the sample molecules in 
the negative direction as separate positive payoffs. This would give the absolute molec- 
ular flow in the positive and negative directions, which can then be subtracted to give 
the net molecular flow. And the zero-variance bias function would then be applicable. 


4J 2 EE 

n=0 

OO 

■E 

n=0 

Hence, 


Minimum-Variance Bias Function 

We can also minimize the variance by using a straightforward calculus of variations 
procedure as follows (see ref. 4). To minimize & as given by equation (58) with 
respect to I n (X n ) subject to the boundary condition of equation (52) and I >0, we ob- 
tain 
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E n< X n> 

V X n> 


oo 

Z D E (X )I (X ) 

n rr n' n v n' 


n=0 


= 0 


where are constants to be determined. This gives 


-[2P(X n )W(X n ) - P 2 (X n )] 


E„(XJ 

Tf- 2 - + D n E n< X n> = 0 

^< X n> 


This can be rewritten as 


- [2P(X n )W(X n ) - P 2 (X„)] 1/Z 

multiplying by E fi (X n ) and integrating. Then using equation (52) to solve for D n , we 
can write the minimized-variance bias function as 


V X n> - 


[2P(X n )W(X n> - P 2 (X n )] 


1/2 


/ E „< X n>[ 2!> < X n> W ( X „> - P 2 (X n )] ^ dX n 


(61) 


Since I n > 0 and real, then [2P(X n )W(X n ) - P 2 (X n )] > 0. If all the P(X n ) > 0, then 
using equation (18) we find 


2P(X n )W(X n )-P 2 (Xj = 


-/l 


n+i |X n )dX n4 , >0 


as required. 

Substitution of I n ( x n ) into equation (56) gives for the minimized expected squared 
bias payoff 



We can show that 



is less than 



as follows: 


Consider the inequality 
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J* ^2P(Y)W( Y) - P 2 (Y)]^ 2 - J |2P(X)W(X) - P^X)^ 2 E n (X)dx| E n (Y)dY > 0 


Then we can show that 



E n (Y)dY 


-VlA 


LuV L 

n=0 


2P(Y)W(Y) 



E n (Y)dY 


r 

}' 


> 0 


This can be written as 


‘MM**) 

This equation shows the Monte Carlo error would be larger for the unbiased as compared 
to the minimized-variance bias case. Since from equation (27) 

2P(X n )W(X n ) - P 2 (X n ) = Q(X n ) -f Q(X n+1 )K(X n+1 |X n )dX n 

The result for the minimized-variance bias function equation (61) implies that the prob- 
ability distributions should be increased at those values of X fi that contribute most to 
the variance. 

a 

The minimized-variance bias function I is not the same as the zero-variance 
t n 

bias function P and so must be a relative minimum because it does not give a zero 
variance as 1^ does. This is true for the case where the Monte Carlo calculation is 
carried out with bias functions exactly as given and where an infinite number of terms 
are evaluated. In the actual calculation, since the bias functions are approximated and 
only a finite number of transitions are evaluated, it is not known which bias function is 
preferable. 


Molecular Diffusion Numerical Results 

We wish to apply the results of the random walk theory to the specific problem of 
rarefied-gas diffusion where the diffusing molecule is considered to be undergoing a 
Markov random walk. The diffusion molecule originates at the emitting plane surface 
with the initial probability distribution given by Eq(Xq). Because the geometry is one 
dimensional (see fig. 1) where the direction from the emitting plate to the parallel ab- 
sorbing wall is along the Z-axis, we need only be concerned with the Z -component of the 
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molecular position. Inasmuch as in this report we will be finding only the molecular flux 
and molecular density at various positions across the gap, we also need consider only 
the Z -component of velocity. The position of the molecule is designated by the dimen- 
sionless coordinate z, which represents the distance that the molecule is from the 
emitting plate divided by the mean free path A . Similarly, the symbol v is employed 
to represent the velocity, where v equals the Z -component of velocity V divided by 
the thermal velocity c. Thus, the state of the molecule is given by the ordered pair 
{z,v}. If we assume thermal equilibrium for the molecules leaving the emitting wall, 

Zq = 0, the initial distribution of the state of the molecule is given by (see appendix A, 
eq- (A2)) 


E 0 ( Z 0 ’V = [ 2v o ex p(- v o)] 6 ( z o> (62) 

The transition probability density can be written as 

K < z n+1" v n+l I z n’ v n> * T(z n+1 1 z n- v n )C < v m-l I z n+l’ V < 63) 

The transport distribution T(z n+ ^ |z , v ) is the probability of having a collision in 
dz , at z , after leaving z . The collision distribution C(v . I z ,, v ) is the 
probability that a molecule with velocity v fi after undergoing a collision at z r+ j will 
have a velocity after collision in dv n+ j at v n+ i- However, once the molecule is inci- 
dent on a wall either at z = 0 or z = l (where l is the inverse Knudsen number, i. e. , 
the distance between the plates divided by the mean free path L/A c ), the molecule his- 
tory is ended. The molecule is then assumed to remain at the wall with zero velocity so 
that the transition kernel is then given by 

K(z n + 1’ v n+l I z n’ V = 6(z n + l " z r) 6 K + J for z n = 0 or 1 ’ n > 0 ^ 

where 6 refers to the usual Dirac delta function. 

As shown in appendix A (eq. (A10)), the transport probability for the case of 
Maxwellian-type collisions, where the collision rate 0 is a constant, can be written as 

T ( z n+1 1 z n> v n )dz n + l = |7-| ex P f | “Vl (65) 

In the present case, the molecules are assumed to emerge from the collision with a 
Maxwellian distribution. Thus, the collision probability is not a function of the previous 
velocity and position and can be written as (see appendix A, eq. (A17)) 
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( 66 ) 


C(v)dv = — — exp(-v^)dv 

(Other types of collisions, either classical or quantum mechanical can be used. ) 

We are interested in determining the macroscopic flow quantities: the molecular 

flux and molecular density of the diffusing species at various positions z across the 

s 

gap between the walls. These positions z are called scoring positions and are illus- 

s 

trated in figure 1. Thus, z = (0. 1)£ corresponds to a scoring position at 1/10 of the 

s 

distance across the gap from the emitting surface. 


Scoring Payoff 

Each sample molecule that crosses the scoring position contributes, to the measure- 
ment of some macroscopic flow quantity of interest, and amount given by an appropriate 
scoring function n. This scoring function will, in general, be a function of the dimen- 
sionless velocity of the molecule as it crosses the scoring position. Since the number 
of molecules per unit time crossing the scoring position in the steady-state case ju 0 is 
directly proportional to the rate at which the molecules leave the emitting plate, pq + , 
the net flux of n crossing the scoring position can be expressed as " 

4 70 _ ^/s[ V7r < v) ] 

M 0+ ^0+^0+^ 


The 2^ n^v) is obtained by scoring the value of n for each sample molecule passing 
i 

the scoring position and then averaging this result over N samples. The i> is the 

local molecular density at the scoring position s, while <^ s E V7r ( v )]] is the value of 

vn(v) averaged over the local molecular velocity distribution of the molecules at the 

scoring position s and is therefore equal to ijS n \ the net mass flux of n transported 

s 

across the scoring position. Similarly, is the molecular density of the molecules 

being emitted from the surface, while <f q + (v) is the value of v averaged over the 
molecular velocity distribution of the molecules at the surface leaving the emitting walls 
per unit time and therefore equals Mq + , the mass flux emitted from the wall. 

When the molecular flux p at the scoring position is desired, the appropriate 
scoring payoff is 


£ *i( v ) (67) 


M 


= ±1 


when 


v 


( 68 ) 
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Then equation (67) becomes 


M _ ^s (v) _ 

^0+ M 0+ M 0+ 



(69) 


When the desired macroscopic quantity is the molecular density v at the scoring posi- 
tion, the appropriate scoring payoff is 



For this case, equation (67) becomes 


(ff ) 

^s V _ ^s^ 1 ) 
^ 0 + ^ 0 + 




(70) 


(71) 


Analog Calculation 

For the first example the "analog" Monte Carlo method is used. First, the sample 
molecule birth velocity Vq is randomly picked using equation (A4). Then similarly a 
point Zj of first collision is obtained from equation (A12). If the sample molecule in 
going from Zq to z^ passes z g , the appropriate scoring payoff 7 t(vq) is scored as 
the zero event payoff Pq. If the sample molecule does not pass z g , = 0. Then a 
new velocity v^ after first collision is found from equation (A21), and a position of 
second collision is found as before. Again if z g is passed by the sample molecule 
in going from Zj to z then P^ = rr(v^); if not, Pj = 0. The process is continued in 
this manner until the molecule is incident on either wall where the random walk termi- 
nates. The sum of all the event payoffs for this random walk then gives t?q. The ran- 
dom walk process is repeated until N independent samples are obtained yielding the 
average value t)q, which is an unbiased and consistent estimate of the expected value of 
the random walk payoff X. 

In figure 2 the molecular flux across the gap per unit molecular flux leaving the sur- 
face is plotted as a function of inverse Knudsen number l . 

In figure 3 is shown the molecular density at the various scoring positions across 
the gap per unit molecular density of the molecular flux entering the gap from the emit- 
ter. These results are shown for a different inverse Knudsen number l . 


22 



In table I the value of the normalized molecular density i' s /i'o+ a * the scoring po- 
sition adjacent to the absorbing wall z = l is given for the analog calculation. Also 

b 

given is the sample standard deviation S(t?q) obtained by use of equation (21). And in 
addition, the average running time per sample a is shown. 


Next -Event Calculation 

In this calculation we take the first step to reduce the variance. As shown previ- 
ously, if the variance can be decreased for each event payoff P , the variance of the 
random walk payoff tiq would also be reduced. In the analog calculation, we registered 
a nonzero event payoff only if a sample molecule actually crossed the scoring position. 

In the present case, called the "next event" calculation, we register a nonzero event 
payoff at each molecular collision within the gap whether or not the sample molecule 

actually crosses the scoring position. We can achieve this as follows: We can write 
th 

the n n expected event payoff as 


X n - fj *<'g E 0< X 0> K<n) < X n |X 0 >-rC Zs , x n ) dx o < 72 > 

where the term r(z , X ) is the scoring probability, that is, the probability that a mole- 
b n 

cule in state X n will reach the scoring position z g without incurring a molecular col- 
lision. Therefore, immediately after each collision we can score the event payoff 

p(x n ) = 7r ( v n ) T(z s’ X n } < 73) 

and obtain an unbiased estimate of X . The scoring probability as shown in appendix A 
(eq. (A13)) is given by 


T (z s , X n ) = exp — j ( 74 > 

The results of using the next-event calculation to find the normalized molecular 
density at the scoring position adjacent to the absorbing wall is shown in table I and cam 
be compared to the results for the analog calculation. The next-event calculations gen- 
erally resulted in a decrease in variance with some increase in running time. This did 
not result in a net saving in computing time as shown by the tabulated values of I, which 
is the percent decrease (+) or increase (-) in computing time for. the biased case com- 
pared to the analog case for the same confidence limits (eq. (17)). The negative values 
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of I indicate an increase in computing time is necessary to obtain the same confidence 
limits in the biased case as in the unbiased case. Careful design of the computer pro- 
gram, which can reduce computing time per sample, would be very effective in reducing 
total computer time needed since large numbers of samples are calculated. 


Birth Biasing 

It is possible to use the zero biasing function given in the analysis by equation (58) 
for the initial event, 


T t p (V 

l 0 -~ 

A 0 

because we can evaluate numerically the equation 


(75) 


E 0 (X 0 )P(X 0 )dX 0 (76) 

For this case, the initial or ' 'birth' ' payoff would be given by P^ = Pq/Iq = *-Q’ thus, we 
can take Aq as our initial payoff irrespective of the initial state Xq found for the 
sample molecule. 

The results for the case where we treat the initial event by the birth bias and use 
the next- event calculation for all subsequent events are shown in table I. For the same 
conditions used in the previous examples in table I, it can be seen from the table that 
there is a decrease in variance for the smaller values of l. There is also a decrease 
in running time. This generally resulted in a savings in computer time, as shown by the 
positive values of I in table I. 


Survival Biasing 

The variance in the random walk payoff of the molecule rj 0 is in part caused by 
the variation in the number of payoffs that are scored before the walk terminates. That 
is, if a sample molecule is incident on one of the walls when the number of previous 
scorings is small, it would be expected that the random walk payoff tjq for this mole- 
cule would be small compared to a sample molecule that had a very large number of 
previous scorings. Hence, if we reduce the variation in the number of collisions, we 
could expect to reduce the variation in t)q. One method of accomplishing the aforesaid 
purpose is to bias the transition probability so that the sample molecule is not permitted 
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to reach an absorbing wall. We can write the survival transport distribution for 
the molecules so that the next collision occurs before the molecule is incident on the 
absorbing surface as 


n(s) ( „ j T ( z nlVl) 

• ^ z nl*n-l' ~'ZT\ 

T S (X n -1.) 


for z n < if 


(77) 


where S£ is either 0 or i 
is less than or greater than 
The normalizing factor 

^ <S)(X n-l> = ^JVl^n <™> 

z n -z n-l 

These results are shown in appendix A (eq. (A15)). 

Since the event payoff P(X ) is zero for all z > if, we can write the expected value 
th ^ 

for the n event payoff as 

A n =/. .y'fpCX^T^^X^T^^Xp . . . T^X^)] 

XE 0 (X 0 )T (s) ( Z;l |x 0 )C( Vl )T (s) (z 2 IX^CCvg) . . . T (s) (z n Ix^C^)] dX Q . . . (79) 

The Monte Carlo procedure for the survival biasing case now consists of picking the new 
position X n by using the survival-biased transport probability T^(z n |X n _^) and scor- 
ing the survival -biased payoff 

P i S) = P(X n )T (s) (X 0 )T (s) (X 1 ) . . . T (s) (X n _ x ) (80) 


depending on whether the sample molecule velocity v n 
0. 

T^ is given by 


Velocity Biasing 

Much of the event payoff variance is caused by the large variations in velocity of 
the scoring sample molecules. In the zero biasing analysis we wished to minimize the 
variation in the event payoff. We can do this by defining a biased collision distribution 

c(C) ( v J Z J = C(v n )l(C)(v n’ z n ) (81) 
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( 82 ) 


where I^ c \x n ) is the collision bias function and is subject to the condition 

/CIVMV z„)dv n - 1 

The biased random walk distribution is now given by 

H (c) = E 0 (X 0 )K (S) (X 0 |X 1 )K^ S )(X 2 |X 1 ) . . . T^ s) ( Zn |X n . 1 )C( Vn )I^ c \x n ) (83) 

and the biased random walk payoff is now 


p( X ) 

t^^XqJt^^Xj) . . . T^ s \x n _ 1 ) 

i (c) (x n ) 


(84) 


The variance of the event payoff °vP n ) is given by 


E ( P n) = -*l-f [pi C) ] 2 H (c) <K 0 ' ' • 


(85) 


We can minimize equation (85) with respect to I^ c \x n ) by using the calculus of variations 
to obtain 



o 

where D is the Lagrange multiplier. This result can be rewritten as 


DI (c) (X n ) = P(X n )T (s) (X 0 )T (s) (X 1 ) . . . T (s) (X n _ x ) (87) 

which can be solved for D by using equation (82) to give 


26 



( 88 ) 


n CP(z n ) 


where 


CPU n ) = I C(v n )P(X n )dv n 
Then the velocity-biased n in event payoff becomes 

P^ C) = PC(z n )T (s) (X 0 )T (s) (X 1 ) . . . T (s) (X n _ x ) (89) 

While the random walk needed for scoring is given by 


H (c) = E 0 (X 0 )T (s) ( Zl |X 0 )C( Vl ) . . . T (s) (z n |X n _ 1 ) (90) 

We can now proceed to use the technique of survival and velocity biasing in evaluating 
Xy Of course, can be evaluated by numerical integration and this value used for all 
first-event payoffs, as in the birth case, regardless of the state X^. However, this nu- 
merical integration would be difficult and instead we proceed as follows: We wish to 
evaluate the biased payoff previously derived in equation (89) as 

P^) =T (s) (X 0 )CP( Zl ) (91) 

The biased random walk is thus given by equation (90) as 

H (c) = E 0 (X 0 )T (s) ( Zl |x 0 ) (92) 

(c) 

To evaluate P^ , we randomly choose Vq from Eq(Xq) as shown in appendix A 

(eq. (A3)). This allows us to evaluate T^ s \xg). Then z^ is found by randomly pick- 
ing T^ s \ Zl |Xq) as given in appendix A (eq. (A16)). With Vq and z^ evaluated, the 
first-event payoff P ^ is scored. We can continue and evaluate P(X2), P(Xg), . . . 
with an unbiased procedure so that the random walk is then given by 

H (c) = E 0 (X 0 )T (s) (z 1 |X 0 )C(v 1 )K(X 2 |x 1 ) . . . (93) 
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Then for the correct 


th 


n 


payoff we have 


P i C) = P( x n )T (s) (X 0 ) n > 1 (94) 

The result for the birth bias and first-term bias with the remaining walk unbiased is 
shown in table I under the ”+1 Term bias. " The results show a significant decrease in 
sample variance, except for the largest value of l, although there is a corresponding 
increase in computer running time. The net savings in computer time was mixed, im- 
proving some of the cases only. 

The biasing can be continued to the second collision, which would have the payoff 


pUO = t ( s )(x 0 )t (s) (x 1 )cp(z 2 ) 


(95) 


The biased random walk is now given by 


H (c) - E 0 (X q )T (s) ( Zi |X 0 )C( Vl )T (s) (z 2 |X x ) 


(96) 


we have previously obtained values of X Q , z ^ and T^ s \v Q ) used in 
finding P^. We then randomly pick v^ from C(vj) and evaluate T^Xj). Then 

>p(s)f 


(z 0 |X«) as shown in appendix A (eq. (A15)). We 
Z ,(c) 


To evaluate 

P S 

is found by randomly picking from 

can then evaluate the biased second-event payoff P 2 W . We can continue in either a 
biased or unbiased manner in this fashion. The result for biasing to the third collision 
P^(Xq) + P^(X^) + P^(X 2 ) + P^ c ^(Xg) while the remaining collisions are unbiased 

T^ s \x 0 )T^ s \x 1 )T^ s \x 2 )[p(X 4 ) + P(X 5 ) + . . .] are shown in table I labeled ”+3 Term 
bias. ,T Now significant reduction in variance is seen even for l of 50 in this case. 
There was a net savings in computer time for all values of l . 


Russian Roulette 

Finally we can bias continuously; however, in this case, the random walk does not 
end because the transport probability biasing does not allow the molecules to be incident 
on the walls. In this case, Russian roulette was used to end the random walks. The 
sample histories were followed until the weighting T^ s ^(Xq)T^ s ^(X 4 ) . . . T^ s \x n _ 4 ) 
was less than 0. 001. Then if a randomly picked number uniform between zero and unity 
was greater than 0. 1, the history was ended. If less, the sample weight was multiplied 
by 10 and the process continued. These results are given in table I under "Russian 
roulette. " This resulted in decreasing the variance but increasing the computer time 
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such that the net saving in time was not greatly improved over the three-term bias case. 
However, no effort was made to optimize the values used, 0. 001 or 0. 1. Further effort 
in this direction could significantly improve the Russian roulette procedure. 


CONCLUSIONS 

Molecular diffusion through a rarefied gas was analyzed by the theory of Markov 
random walks. The results indicate variance reduction techniques can be successfully 
used in Monte Carlo analyses of rarefied-gas problems. This can be important in cases 
where "analog" Monte Carlo analysis of unbiased random walks gives such large vari- 
ances that unreasonably large amounts of computer time are needed to obtain acceptable 
confidence limits on the results. 

In the present analysis, both a zero-variance and a minimum-variance biasing func- 
tion are found. These biasing functions, however, apply only for positive event payoffs. 
However, by appropriately defining the scoring, the event payoffs can be made positive. 

The minimum -variance biasing function does not give a zero variance, as the zero- 
variance biasing function does, and so must give a relative minimum variance. How- 
ever, in the Monte Carlo calculation procedure when only a finite number of events can 
be calculated and the bias functions are approximated, which of the biasing functions is 
optimum has not been determined. 

The present results can be applied to any Markov process such as thermal radia- 
tion, where instead of following molecules through various molecular collisions we can 
follow photons through the random walks and apply these same variance reduction tech- 
niques. 

The "Russian roulette" procedure did not significantly increase the running time 
saved. The reason may be that the parameters used in the Russian roulette procedures 
were not optimized. 

Since the Monte Carlo procedure is a continuous sampling process, further improve- 
ment can be obtained in computer time saving by introducing a learning process into the 
program so that the bias function would be changed as additional samples are obtained. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, June 22, 1973, 

502-28. 
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APPENDIX A 


BIRTH DISTRIBUTION 

Assume that the molecules leaving the emitting surface at z = 0 are in thermal 
equilibrium with the emitting surface and that the z -component of velocity v in dv is 
divided by the total molecular flux leaving the surface jliq + given by 



(Al) 


Then the flux of molecules leaving the emitting surface at z = 0 gives the initial 
z-component of velocity distribution of the molecules leaving the surface (ref. 3) as 

Eq(v) = 2v exp(-v 2 ) (A2) 

In the Monte Carlo calculation we can randomly choose a Vq from this distribution 
by using 


2v z exp(-v z jdv z = 1 - R = R (A3) 

where R is a random number generated on the computer from a uniform distribution 
between zero and unity. Equation (A3) can be rewritten as 

v Q = (-In R) 1/2 (A4) 

Path Length to Collision 

The distance a molecule travels between collisions can be found as follows: The 
probability of a molecule having a collision in time dt is given by 

- = +© dt (A5) 

N 

where 0 is the collision rate. The probability of not having a collision during time t 
is given by N/Nq, which can be obtained by integration of equation (A5) to give 
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N 


N, 


= exp(-0t) 


0 


(A6) 


Then the probability for having a collision in dt after having no collision in time t is 
given by 


[exp(-0t c )j (-©)dt 


The distance between collisions is then given by 

A = |v|t c =M ( t c 0) 


(A7) 


(A8) 


In reference 3 for a Maxwellian-type collision the collision rate 0 is shown to be a 
constant. Then if we define a mean free path X £ = C/0 and normalize so that 6 = A/X c 
we obtain for the distance between collisions 


6 =A=M( t 0) = | v | (t 0) (A9) 

X c C 

Combining equations (A7) and (A9) gives the probability of having a collision to be 


T<z n+ll z n v n ) = 


exp - 


n 


n' 


d6_ 


n 1 


where 


6 = z i - z 
n n+1 n 


We can randomly pick from this distribution by using 

• 6 . 


R 


= / " T(6|v s 


)d6 


This can be integrated to give 


6 n = -|v z |lnR 


(A10) 


(All) 


(A12) 
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The scoring probability, which is the probability of reaching a scoring position z 

b 

from some point z without having an intervening collision, is given by the probability 
of having a first collision at 6 > 6 = I z - z | . This is given by 

b & 


t( 6 s | v z ) = 


r 

y 6=6. 


T(6 |v)d6 = e 


■ 6 ./ 1 1 


(A13) 


In some biasing cases we wish to pick from a transport distribution containing dis- 
tances to collision which are limited to values smaller than those which would cause 
the molecules to be absorbed at the walls: 


5 ^ where 



z | for v > 0 
for v < 0 


For this case the transport probability is given by 



We can randomly pickfrom this distribution by using 


6 




(A14) 


(A15) 


(A16) 


Sample Velocity After Collision 

The molecules coming out of collision are assumed to be in Maxwellian equilibrium. 
We can thus write for the distribution of molecules coming out of collision the Maxwel- 
lian v-component of velocity as 


1 " 

C(v)dv = — — e v dv 
Vtt 


(A17) 


32 



To randomly sample from this distribution we can write 


2 -X^ 2 

f(v, X) = e” v — dvdX = 2e' p p dp — dp dd 

y 77 yfn n 


We can then find 


so that 


v = p cos 9 0 < p < 00 

X = p sin 0 0 < 0 < 7r 

p = (-In R^ 1 / 2 

0 = Rgff 

v = (-In Rj) 1//2 cos(^R 2 ) 


(A19) 


(A20) 


(A21) 
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APPENDIX B 


SYMBOLS 


c <VllVl’ v n> 

c(c) (v n |z n ) 

C 

D 

E(X) 

E n (X) 

E 0 (X) 

4 ) 

' 8 (> 

g(x|y) 

H(X 0 ,X 1S . . .) 

I 

x (c) 


K(Y|X) 

K^)(Y|x) 

<£ 

L 

L n 

l 

P 

Q(X) 

R 


collision distribution 

biased collision distribution 
thermal velocity 
constant 

distribution at X 
distribution at X after n events 
initial or birth distribution 
expected value 

expected value based on local molecular distribution 

Green's function 

random walk distribution 

percent running time (saved +, lost -) 

collision bias function 

bias function for n events 

zero -variance bias function 

minimized-variance bias function 

transition distribution 

transition distribution for l events 

limiting value 

distance across gap 

biasing ratio 

inverse Knudsen number, L>A C 
event payoff 

conditional expected square payoff for a random walk leaving X 

random number picked from a uniform distribution between zero 
and unity 

sample variance 
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T ‘ z nl x n-1> 

T (S, (Z n l X n-l> 

T <S) (X n -i> 

V 
v 

W(X) 

X n 

Y 

Z 


a 

A 

6 


e 


^0 

0 

X 


X A 
ju’ n 


M 




(*r) 


^0+ 

n 


a 


T 

Subscripts: 

c 


transport distribution 
survival transport distribution 
normalizing factor 

point in velocity space (Z -component of velocity) 
dimensionless velocity in z-direction, V/c 
conditional expected payoff for a random walk leaving X 

1L 

point in position and velocity space immediately after n n 

point in position space 

coordinate across gap 

dimensionless distance across gap, Z/X c 

average running time per sample 

distance between collisions 

dimensionless distance between collisions 

confidence limit 

random walk payoff 

collision rate 

expected payoff for random walk 
mean free path 

expected value of k th event payoff 
expected molecular flux and density payoff 
molecular flux 

molecular flux of it at scoring position s 
molecular flux emitted at plate 
molecular density at scoring position 
scoring function 
variance 

scoring probability 
thermal velocity 


event 
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n event index 

s scoring position 

z position in channel 

Id molecular flux 

u molecular density 

0 evaluated at z = 0 

+, - positive and negative z- direction 

Superscripts: 

(s), (c) survival biasing, velocity biasing 
Id molecular flux 

v molecular density 

t zero-variance case 

' minimum -variance case 

* biased case 

sample average 
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TABLE I. - NUMERICAL RESULTS - 10 000 SAMPLES 


Monte Carlo 
method 

Density 

ratio, 

u/v 0+ 

Standard 

deviation, 

S(no) 

Average 
running 
time per 
sample, 

a, 

min 

Percent 
change in 
running 
time 
(saved +, 
lost -), 

I 


l 

= L/A c = 0. 1 



Analog 

0.848 

0. 923 

0. 13X10' 4 

0 

Next event 

. 844 

. 867 

. 15 

-1.8 

Birth bias 

. 845 

.741 

. 12 

+40. 5 

+1 Term bias 

. 850 

. 1568 

.42 

+90.6 

+3 Term bias 

. 855 

. 116 

.64 

+92.2 

Russian roulette 

. 854 

. 115 

.61 

+92.7 


l 

= L A C = 1 



Analog 

0.4914 

0.7494 

0. 31X10' 4 

0 

Next event 

.4881 

.7188 

. 35 

-3.8 

Birth bias 

.4881 

.7185 

.32 

+5. 1 

+1 Term bias 

.4958 

. 5728 

. 59 

-11.2 

+3 Term bias 

.487 

.336 

.95 

+38.4 

Russian roulette 

. 491 

.259 

2.33 

+10.2 


l = 

L/A c = 10 



Analog 

0. 1202 

0. 578 

2. 14X10" 4 

0 

Next event 

. 11965 

. 527 

2. 33 

+9. 49 

Birth bias 

. 11965 

. 527 

2. 31 

+ 10.2 

+1 Term bias 

. 1171 

. 463 

2. 48 

+25.6 

+3 Term bias 

. 1195 

. 402 

3. 22 

+27.2 

Russian roulette 

. 1172 

. 244 

13. 97 

-16. 3 


l = 

L/\ c = 50 



Analog 

0.0261 

0. 224 

10. 78X10' 4 

0 

Next event 

.0263 

. 238 

11. 61 

-21. 6 

Birth bias 

. 0263 

. 238 

11.60 

-21. 5 

+1 Term bias 

.0263 

. 2759 

11. 34 

-59. 6 

+3 Term bias 

.0246 

. 183 

14. 17 

+ 12.27 

Russian roulette 

. 0245 

. 168 

39. 50 

-106. 1 




• Nondiffusing species 

O Diffusing species 

Scoring cross sections, Z s 




Figure 2. - Molecular flux across gap per unit molecular 
flux as function of inverse Knudsen number - analog 
calculation for 10 000 samples. 


Inverse Knudsen 



Dimensionless distance across channel, z = Z/L 


Figure 3. - Molecular density at various scoring positions 
across gap per unit molecular density entering gap from 
emitter - analog calculation for 10 000 samples. 
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